Geologists have discovered that one can identify certain physical characteristics of a geological formation from the gravitational potential field (denoted with the symbol G in this application) near the formation. For example, the gravitational field G can often indicate the presence and yield the identity of a mineral such as coal that is located beneath the surface of the formation. Therefore, measuring and analyzing the gravitational field G of a formation, can often yield such physical characteristics of the formation more easily and less expensively than an invasive technique such as drilling. Relevant characteristics of the field are typically determined not by measuring the gravitational potential G directly, but either by measuring components of the gravitational acceleration vector g resulting from this field or by measuring spatial derivatives of these acceleration vector components. The three components of the acceleration vector can each be spatially differentiated along three different axes providing a set of nine different signals mathematically related to the underlying gravitational potential G. These nine signals are the gravitation tensor elements (sometimes called the gravity Gradients), and much effort has gone into developing techniques for accurately measuring these tensor elements.
Referring to FIG. 1, one can use a gravity gradiometer 10 to measure the gravitational potential field G near a geological formation (not shown). The notation used in this patent to the refer to the nine gravitational field tensor elements in matrix form is:                     Γ        =                  [                                                                      Γ                  xx                                                                              Γ                  xy                                                                              Γ                  xz                                                                                                      Γ                  yx                                                                              Γ                  yy                                                                              Γ                  yz                                                                                                      Γ                  zx                                                                              Γ                  zy                                                                              Γ                  zz                                                              ]                                    (        1        )            
where the matrix members represent the respective gravity tensors along each of the three X, Y, and Z xe2x80x9cbodyxe2x80x9d axes, which typically intersect at the centroid 12 of the gradiometer 10. For example, the tensor element xcex93xx (which can be expressed in equivalent units of (meters/seconds2)/meter, 1/seconds2, or Exc3x6tvxc3x6s units where 109 Exc3x6tvxc3x6s=1/seconds2) is the spatial partial derivative along the X axis of the X component of the gravitational acceleration vector g, xcex93xy is the spatial partial derivative along the Y axis of the X component of g, xcex93xz is the spatial partial derivative along the Z axis of the X component of g, xcex93yx is the spatial partial derivative along the X axis of the Y component of g, etc. Furthermore, although the tensors elements xcex93 may vary over time, for many formations the tensor elements xcex93are constant over time, or vary so slowly that they can be treated as being constant over time. Moreover, in some applications the gradiometer 10 may make measurements sufficient to calculate only the desired elements of the full tensor xcex93. To measure the gravitational potential field G of a geological formation (not shown in FIG. 1), one mounts the gradiometer 10 in a vehicle (not shown) such as a helicopter that sweeps the gradiometer over the formation. For maximum accuracy, it is desired that the gradiometer 10 not rotate at a high rate about any of the X, Y, and Z body axes as it sweeps over the formation. But unfortunately, the vehicle often generates vibrations (e.g., engine) or is subject to vibrations (e.g., wind) that causes such rotations about the body axes. Therefore, the gradiometer 10 is often rotationally isolated from the vehicle by a gimballing system (not shown) which allows the gradiometer to remain non-rotating even as the vehicle experiences varying orientations typical of its operation. The gimballing system carrying the gradiometer 10 typically includes a rotational sensor assembly 18, such as a gyroscope assembly for measuring rotational activity (typically rotational rate xcfx89 or rotational displacement) about the X, Y, and Z body axes. Control signals derived form these measurements are fed back to the motors attached to the gimbal axes to reduce the rotations experienced by the gradiometer 10. But although the gimballing system typically reduces the magnitude of the vibration-induced rotations of the gradiometer 10 about the body axes, it is typically impossible to eliminate such rotations altogether. Even with an ideal gradiometer, the tensor measurement would, of physical necessity, be additively corrupted by the presence of gradient signals caused by these rotations. These additional non-gravitational gradients are simple deterministic functions of the rotational rates (e.g. rotational xcex93=xcfx89x, xcfx89y where xcfx89j refers to the rotational rate around the j body axis in radians/sec). Consequently, the measurements from the gradiometer 10 typically have these corrupting signals subtracted by a processor 20 to increase the accuracy of the gradiometer""s measurement of the gravitational field G as discussed below in conjunction with FIG. 3. Although shown as being disposed within the housing 16, the processor 20 may be disposed outside of the housing for processing of the measurement data from the gradiometer 10 in real time or after the gradiometer measures the gravitational field G. In the latter case, the gradiometer 10 typically includes a memory 22 for storage of the measurement data for later download to the external processor 20. Alternatively, the gradiometer 10 may include a transmitter (not shown) for transmitting the measurement data to the external processor 20 and/or an external memory 22. Moreover, the processor 20 or memory 22 includes a sample-and-hold circuit (not shown) and an analog-to-digital converter (ADC) (not shown) to digitize the gradiometer measurement data and any other measured signals required for optimal operation.
Referring to FIG. 2, the gravity gradiometer 10 of FIG. 1 includes one or more disc assembliesxe2x80x94here three disc assemblies 24, 26, and 28xe2x80x94each for measuring a subset of the full set of tensors F for the gravitational field G of a geological formation 36.
Each disc assembly 24, 26, and 28 includes a respective disc 30, 32, and 34 that is mounted in a respective plane that is coincident or parallel with one of the three body-axis planes such that the spin axis of the disc is either coincident with or parallel to the body axis that is normal to the mounting plane. Furthermore, each disc includes orthogonal disc axes that lie in but rotate relative to the mounting plane. For example, the disc 30 lies in the X-Y body-axis plane, has a spin axis Zs that is parallel to the Z body axisxe2x80x94that is, the X-Y coordinates of ZS are (X=C1, Y=C2) where C1 and C2 are constantsxe2x80x94and includes orthogonal disc axes XD and YD. As the disc 30 rotatesxe2x80x94here in a counterclockwise directionxe2x80x94the XD and YD, disc axes rotate relative to the non-rotating X and Y body axes. At the instant of time represented in FIG. 2, the XD and YD disc axes of the disc 30 are respectively parallel and coincident with the X and Y body axes. In addition, the disc 32 lies in a plane that is parallel to the Y-Z body-axis plane and has a spin axis XS that is parallel to the X body axis. At the instant of time represented in FIG. 2, the YD and ZD disc axes of the disc 32 are respectively parallel with the Y and Z body axes.
To measure the gravitational field G, the disc assemblies 24, 26, and 28 each include at least one respective pair of accelerometers that are mounted xcfx80 radians apart on the discs 30, 32, and 34, respectively. For clarity of explanation, only the disc assembly 24 is discussed, it being understood that the other disc assemblies 26 and 28 are similar. Here, the disc assembly 28 includes two pairs of accelerometers 38a, 38b and 38c, 38d. Each accelerometer 38a, 38b, 38c, and 38d includes a respective input axis 40a, 40b, 40c, and 40d along which the accelerometer measures a respective acceleration magnitude Aa, Ab, Ac, and Ad, and each accelerometer is mounted to the disc 30 such that its input axis is a radius R from the spin axis ZS and is perpendicular to R. The accelerometers 38a and 38b of the first pair are mounted xcfx80 radians apart on the XD disc axis, and the accelerometers 38c and 38d are mounted xcfx80 radians apart on the YD disc axis. Although ideally described as being perpendicular to the radius R, the input axes 40a, 40b, 40c, and 40d may actually be oriented at other angles relative to R intentionally or due to manufacturing imperfections. Furthermore, the disc assembly 24 may include additional pairs of accelerometers that are mounted on the disc 30 between the accelerometers 38a, 38b, 38c, and 38d. For example, the disc assembly 24 may include additional accelerometers 38e, 38f, 38g, and 38h, which are respectively spaced xcfx80/4 radians from the accelerometers 38a, 38b, 38c, and 38d. As is known, these additional accelerometers allow, through redundant measurement, an increase in the signal-to-noise ratio (SNR) of the gravitational-field measurement.
Referring to FIG. 3, the operation of the disc assembly 24 is discussed, it being understood that the operation of the disc assemblies 26 and 28 of FIG. 2 is similar.
FIG. 3 is a top view of the disc assembly 24 where the spin axis Zs extends out of the page from the center 50 of the disc 30. For purposes of explanation, the following ideal conditions are assumed. First, the disc 30 spins in a counterclockwise direction at a constant rate of xcexa9, which has units of radians/second. Second, the input axes 40 are each aligned perfectly with either the XD or YD disc axes and as a result lie in or parallel to the X-Y plane. Third, the accelerometers are all the same radial distance R from the ZS spin axis. And fourth, there are no rotations of the disc 30 about the X or Y body axes.
At a time t=0, the XD and YD disc axes of the disc 30 are respectively parallel and coincident with the body axes X and Y. As the disc 30 spins, the XD disc axis forms an angle xcexa9t relative to its initial (t=0) position. To illustrate this rotation, the position of the XD and YD axes and the accelerometer 38a are shown in dashed line at xcexa9t=xcfx80/4 radians. Although not shown in dashed line, the other accelerometers 38b, 38c, and 38d are also xcfx80/4 radians from their illustrated (xcexa9t=0) positions when xcexa9t=xcfx80/4 radians. Consequently, an equation that represents the acceleration Aa in terms of the gravity tensor elements xcex93xx, xcex93xy, xcex93yx, and xcex93yy, can be derived as follows, where ax and ay are the gravitational-field induced accelerations at the center 50 in the X and Y directions, respectively. Specifically, Aa equals the component of acceleration along the input axis 40a caused by an acceleration in the Y direction minus the component of acceleration along the input axis caused by an acceleration in the X direction. Therefore,
Aa=(ay+xcex93yxR cos xcexa9t+xcex93yyR sin xcexa9t) cos xcexa9txe2x88x92(ax+xcex93xxR cos xcexa9t+xcex93xyR sin xcexa9t) sin xcexa9txe2x80x83xe2x80x83(2)
Expanding the terms of equation (2) gives:
Aa=ay cos xcexa9t+xcex93yxR cos2t+xcex93yyR sin xcexa9t cos xcexa9txe2x88x92ax sin xcexa9txe2x88x92xcex93xxR cos xcexa9t sin xcexa9txe2x88x92xcex93xyR sin2xcexa9txe2x80x83xe2x80x83(3)
Using trigonometric identities for cos2xcexa9t and cos xcexa9t sin xcexa9t, and realizing that xcex93xy=xcex93yx for any gravitational field G, one obtains:                     Aa        =                                            a              y                        ⁢            cos            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    -                                    a              x                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    +                                    Γ              xy                        ⁢                          R              ⁡                              (                                                      1                    2                                    +                                                            1                      2                                        ⁢                    cos                    ⁢                                          xe2x80x83                                        ⁢                    2                    ⁢                    Ω                    ⁢                                          xe2x80x83                                        ⁢                    t                                                  )                                              +                                    Γ              yy                        ⁢                          R              2                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    -                                    Γ              xx                        ⁢                          R              2                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    -                                    Γ              xy                        ⁢                          R              ⁡                              (                                                      1                    2                                    -                                                            1                      2                                        ⁢                    cos                    ⁢                                          xe2x80x83                                        ⁢                    2                    ⁢                    Ω                    ⁢                                          xe2x80x83                                        ⁢                    t                                                  )                                                                        (        4        )            
And combining terms of equation (4) gives:                     Aa        =                                            a              y                        ⁢            cos            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    -                                    a              x                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    +                                    Γ              xy                        ⁢            R            ⁢                          xe2x80x83                        ⁢            cos            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    +                                    R              2                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢                          t              ⁡                              (                                                      Γ                    yy                                    -                                      Γ                    xx                                                  )                                                                        (        5        )            
Because the accelerometer 38b is always xcfx80 radians from the accelerometer 38a, one can easily derive the following equation for the acceleration magnitude Ab by replacing xe2x80x9cxcexa9txe2x80x9d with xe2x80x9cxcexa9t+xcfx80xe2x80x9d in equations (2)-(5):                     Ab        =                                            -                              a                y                                      ⁢            cos            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    +                                    a              x                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    +                      R            ⁢                          xe2x80x83                        ⁢            cos            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t            ⁢                          xe2x80x83                        ⁢                          Γ              xy                                +                                    R              2                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢                          t              ⁡                              (                                                      Γ                    yy                                    -                                      Γ                    xx                                                  )                                                                        (        6        )            
Summing equations (5) and (6) gives the following equation for the expected idealized output of the sum of these two accelerometers:
Aa+Ab=2xcex93xyR cos 2xcexa9t+R sin 2xcexa9t(xcex93yyxe2x88x92xcex93xx)xe2x80x83xe2x80x83(7)
To increase the accuracy of the measurement (in view of potential errors as discussed below), one can derive equations that represent Ac and Ad in terms of the gravity tensor elements xcex93xx, xcex93xy=xcex93yx, and xcex93yy by respectively replacing xe2x80x9cxcexa9txe2x80x9d with xe2x80x9cxcexa9t+xcfx80/2xe2x80x9d and xe2x80x9cxcexa9t+3xcfx80/2xe2x80x9d in equations (2)-(6) to arrive at the following equation:
Ac+Ad=xe2x88x922xcex93xyR cos 2xcexa9txe2x88x92R sin 2xcexa9t(xcex93yyxe2x88x92xcex93xx)xe2x80x83xe2x80x83(8)
Subtracting these two accelerometer sums (reflected in equations 7 and 8) provides the following equation, which is the basic element of measurement for gradiometers of this design:                                                         (                              Aa                +                Ab                            )                        -                          (                              Ac                +                Ad                            )                                2                =                              2            ⁢                          Γ              xy                        ⁢            R            ⁢                          xe2x80x83                        ⁢            cos            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢            t                    ⁢                      xe2x80x83                    +                      R            ⁢                          xe2x80x83                        ⁢            sin            ⁢                          xe2x80x83                        ⁢            2            ⁢            Ω            ⁢                          xe2x80x83                        ⁢                          t              ⁡                              (                                                      Γ                    yy                                    -                                      Γ                    xx                                                  )                                                                        (        9        )            
This combination signal, which is normally called the bandpass signal, is typically bandpass filtered and digitized, and is then synchronously demodulated by the processor 20 at sin 2xcexa9t and cos 2xcexa9t to recover xcex93xy=xcex93yx and (xcex93yyxe2x88x92xcex93xx)/2.
Still referring to FIG. 3, the conditions that were assumed to be ideal for the derivation of equations (2)-(9) are, unfortunately, seldom ideal. Consequently, these non-ideal conditions introduce additional acceleration terms into these equations, and these terms can reduce the accuracy of the calculated gravity tensor elements if unaccounted for. But fortunately, the processor 20 can account for many of these additional terms as discussed below.
For example, still referring to FIG. 3, the motor (not shown) that spins the disc 30 may be unable to maintain a constant rotation rate xcexa9. Such uneven rotation may cause the pairs of accelerometers to sense reinforcing accelerations that swamp out the accelerations caused by the gravitational field. Consequently, the gradiometer 10 (FIG. 1) may include a sensor (not shown) that measures the rotation rate xcexa9 as a function of time, and the processor 20 can use this measurement to conventionally include the acceleration term introduced by the uneven rotation in the equations (2)-(9).
Furthermore, as discussed above in conjunction with FIG. 1, vibrations of the vehicle (not shown) or other forces may cause the gradiometer 10 to rotate about the X or Y body axes. Such rotations may cause the pairs of accelerometers to sense reinforcing accelerations that swamp out the accelerations caused by the gravitational field. For example, assume that the gradiometer 10 rotates about the Y body axis with a rotational rate (in units of radians/second) xcfx89y. This rotation causes the accelerometer 38a to sense a centripetal acceleration toward the Y axis along a moment arm 52 according to the following equation, where AaY is the acceleration term added to equation (2) due to this centripetal acceleration:
Aay=(xcfx89y)2R cos xcexa9t sin xcexa9txe2x80x83xe2x80x83(10)
The accelerometer 38b senses an identical centripetal acceleration AbY, and the corresponding centripetal accelerations AcY and AdY sensed by the accelerometers 38c and 38d are given by an equation similar to equation (10). Consequently, the processor 20 can use the signals (from the gyroscope assembly 18 of FIG. 1) that represent the rotational rates xcfx89x, xcfx89y, and xcfx89z to include AaY, AbY, AcY, and AdY in equation (9), and thus to compensate the gradient measurements for centripetally induced errors introduced by rotations about the X, Y, and Z body axes.
Similarly, the processor 20 can often account for errors introduced by the input axes 40 of the accelerometer pairs not making the same angle with the respective XD or YD disc axis or not being the same radial distance R from the disc center 50. In these cases, the precise amount of misalignment or radial distance error is typically unknown (though it may be relatively constant for a given gradiometer instrument and may be imperfectly known), and thus the error introduced into the gravitational-field measurement may not be exactly known. However if the functional relationship between the causative error and the resulting signal corruption is known, then this information can be included in an estimation procedure that allows test measurements to be processed, an optimal fit to be made between the measurements identified as being corrupt, and finally a compensation to be applied using these optimal fit estimates. Most often there is a linear (or linearizable) relationship between the error parameters and resultant signal corruptions and a standard least-squares fit is made between the corrupted measurements and arbitrarily scaled calculations of the expected signal corruptions. These expected functions are called regressors, and the fit procedure calculates the extent to which these regressors appear in the raw measurements.
Unfortunately, no set of regressors fits all the acceleration and rotationally induced error found in a gradiometer system. An important part of improving gradiometer instrument performance is the identification of error sources, estimation and compensation of the error effect found in particular instruments, and perhaps adjustment of the instrument build/setup to reduce the physical effects leading to the errors.
Embodiments of the invention as discussed below concern the discovery of one such error mechanism, the calculation of the error effect resulting from this (thus allowing for compensation and improved measurement performance), and identification of instrument adjustments for reducing the magnitude of the raw effect of the error
In one aspect of the invention, a method includes measuring an acceleration along an input axis of an accelerometer mounted to a gradiometer disc, the accelerometer having a coordinate axis that is parallel to a spin axis of the disc, and includes calculating a gravity tensor element as a function of the measured acceleration and a component of the measured acceleration caused by an acceleration along the coordinate axis.
This technique typically yields a more accurate calculation of the gravitational field by accounting for undesired accelerations picked up by accelerometers having input axes that are not parallel to the gradiometer disc. Furthermore, this technique is applicable to systems that measure the full gravitational tensor as well as to those that measure a subset of the full tensor.